function savings_EWM_onedim = generate_results_onedimension(bsimax,sample,wave,tcost,covariate,winter_prevuse,SMC,fullsample,savedir,fbase)
% This function calculates savings from one-dimensional rule based on
% average baseline consumption only. The average baseline consumption
% measure is discretized and searched over each level. 

% Inputs: 
% (1) bsimax: number of bootstrap runs
% (2) sample: cross-sectional dataset merged using elec_wave367_cross.csv
% and propensity scores and renamed as data_estimation_pooled_ps
% (3) wave: indicates whether the output is wave-specific (=3,6,or 7) or
% for the pooled sample (=0)
% (4) tcost: private marginal cost of implementing the program per household
%  >0 represents cost savings; =0 represents kWh reduction 
% (5) covariate: adapted from two-dimensional algorithm, use
% income as the placeholder. 
% (6) winter_prevuse: indicates whether baseline consumption is calculated 
% as the mean of consumption in winter months (Jan and Feb) or as the mean
% of specified pre-treatment periods.
% (7) SMC: =1 use social marginal cost; =0 use retail electricity price
% (8) fullsample: =1 if fullsample, =0 if estimation sample
% (9) savedir: directory to save the savings table 
% (10) fbase: string used to indicate the propensity score and baseline
% months specification for output table 

% Output:
% (1) savings_EWM_onedim: a table summarizing the savings point estimate,
% share of households to treat, total savings for the entire state, and 95%
% confidence intervals for the point estimate 
%% Gather inputs

[X,Y,D,n,ps]=generate_inputs_quadrant_rule(sample,wave,tcost,covariate,SMC);

%% Estimate g using IPW
g = Y.*(D./ps - (1-D)./(1-ps)); % elements of W(G) - W(G_0)

%% Level search
tic1= tic;

% sort rows of X -> [Xr, gr]
[Xr, Ind] = sortrows(X);
gr = g(Ind);

% descretize baseline usage
dx1=10;
[N,boo1,xx1]=histcounts(Xr(:,2),...
    'BinWidth',[dx1],'Normalization','count'); %

nu = length(boo1);

Xu=nan(nu,1); gu=nan(nu,1); nw = nan(nu,1);
k=1;
for i = 1:length(boo1)
       bin_boo = logical((xx1==i));
       % identify individuals in each level of baseline consumption
       Xu(k,:)=boo1(i);
       gu(k,1) = sum(gr(bin_boo)); % net savings for each unique covariate combo
       nw(k,1) = sum(bin_boo); % count number of individuals in the pair of covariate combo
       k=k+1;
end

x1 = Xu(:,1);

% Create grids of midpoints between distinct values of x1/x2
u_x1 = boo1';

eps_x1 = min(u_x1(2:end) - u_x1(1:end-1))/2;
grid_x1 = [u_x1(1)-eps_x1; (u_x1(2:end) + u_x1(1:end-1))/2; u_x1(end)+eps_x1];
k1 = size(grid_x1,1);

% Grid search 
tic1= tic;
m_W = Inf*ones(k1);
v_x1 = grid_x1;

minW = Inf;

for i1=1:k1
    boo_i1=(x1 - grid_x1(i1));
        for sign1 = -1:2:1 %for i = [-1,1]
                select = ((sign1.*boo_i1)>=0); 
          
                %must meet both threshold conditions to be selected as treated
                W = sum(gu.*select);
                
                if W<m_W(i1)
                   m_W(i1) = W; 
                end
                
                if (W < minW)
                    min_i1 = i1;
                    min_sign1 = sign1;
                    minW = W;        
                end
        end
end

fprintf('Serach takes %.2f sec\n',toc(tic1))

%% Save coefficients and other variables for graphs

in_Ghat = (min_sign1*(x1 - grid_x1(min_i1))>=0);

%% Confidence interval
rng(0);

if (tcost)
    if (SMC)
        s_speci = 'smc';
    else
        s_speci = 'pmc';
    end
else
    s_speci='kwh';
end

vhats_p = zeros(bsimax,1); 
vhats_n = zeros(bsimax,1); 
scaled_ATT_itr = zeros(bsimax,1);
ATE_itr= zeros(bsimax,1);
W_bs_Ghat_itr = zeros(bsimax,1);
    
    % shared name for each bootstrap batch
    s_now_batch=datestr(now(),'yyyymmdd_HHMMSS');
    % Create a saving directory for the i_BS files
    fpath_BS = ['.\intermediate_data\EWM_results\BS\quadrant\',...
        s_now_batch,'_','univariate','_',s_speci,'\'];
    if ~exist(fpath_BS,'dir')
       mkdir(fpath_BS); 
    end

tic
    dt_start = now();
for bsi = 1:bsimax
    % DRAW BOOTSTRAP INDEXES
    bsperm = randi(n,[n 1]);
    % Subtract Wn from Wbs
    gr = [g(bsperm,:); -g];
    Xr = [X(bsperm,:); X];
    % Compress the data with identical X's
    [Xr, Ind_bs] = sortrows(Xr);
    gr = gr(Ind_bs);
    [N,coo1,xx1_bs]=histcounts(Xr(:,2),...
    'BinWidth',[dx1],'Normalization','count'); %
    bnu = length(coo1);
    bXu=nan(bnu,1); bgu=nan(bnu,1); bnw = nan(bnu,1);
    kk=1;
    for i = 1:length(coo1)
        bin_coo = logical(xx1_bs==i);
        bXu(kk,:)=coo1(i);
        bgu(kk,1) = sum(gr(bin_coo));
        bnw(kk,1) = sum(bin_coo);
        kk=kk+1;
    end
  
    % Maximize and minimize Wn-Wbs
    maxWb = -Inf; minWb = Inf;
    
    
    for i1=1:k1
        boo_i1=(x1 - grid_x1(i1));
            for sign1 = -1:2:1 %for i = [-1,1]
                    select = ((sign1.*boo_i1)>=0); 
                    Wb = sum(bgu.*select);
                    maxWb = max(maxWb,Wb);
                    minWb = min(minWb,Wb);
            end
    end
    vhats_p(bsi) = maxWb;
    vhats_n(bsi) = minWb;
    sprintf("Iteration %d completed after %0.2g minutes", bsi, toc/60)
    dt_cmpt = now();
        
    %% calculate scaled ATT and ATE for each bootstrap sample 
    scaled_ATT_bs= (sum(D(bsperm,:))/n).*...
        (sum( D(bsperm,:).*Y(bsperm,:) - ((1-D(bsperm,:)).*Y(bsperm,:).*ps(bsperm,:))./(1-ps(bsperm,:)) )/sum(D(bsperm,:)));
    ATE_bs = mean(g(bsperm,:)); % equiv. to ATE = mean(D.*Y./ps)-mean(((1-D).*Y)./(1-ps));
    
    %% apply G_hat to bootstrap sample 
    in_Ghat_on_bs = (min_sign1*(bXu(:,1) - grid_x1(min_i1))>=0);
    W_bs_Ghat = sum(in_Ghat_on_bs.*bgu)/n;
    
    scaled_ATT_itr(bsi,1) = scaled_ATT_bs;
    ATE_itr(bsi,1) = ATE_bs;
    W_bs_Ghat_itr(bsi,1) = W_bs_Ghat;
    
    %% save each BS result
    dt_outp = now();
    dt=[dt_start,dt_cmpt,dt_outp];
    save([fpath_BS,sprintf('%s_BS_%d.mat',s_now_batch,bsi)],...
        'bsi','s_now_batch','Wb','maxWb','minWb','scaled_ATT_bs','ATE_bs','W_bs_Ghat',...
        'bsperm','dt','covariate');
    % progress markers
    fprintf('Boot #%d took %2.3f sec\n',bsi,toc(tic1))
    
end

% two-sided 95% confidence interval for the welfare gain
vhats_abs = max([vhats_p -vhats_n],[],2);

 %% Export rule parameters for graphs 

switch winter_prevuse % Use month_since_opower=[-21,-1] as pre-treatment months 
    case 1
        s_winter = '_winter';
    case 0
        s_winter = '';
end
if (tcost)>0 % cost_savings
    if (SMC)
        s_cost = 'SMC';
    else
        s_cost = 'PMC';
    end
elseif (tcost)==0
    s_cost='kwh';
end

if fullsample==0
    filename_coefs=sprintf('coef_onedim_%s_baseline_%s%s_wave%1.0f.mat',covariate,s_cost,s_winter,wave);
else
    filename_coefs=sprintf('coef_onedim_fullsample_%s_baseline_%s%s_wave%1.0f.mat',covariate,s_cost,s_winter,wave);
end    

filename = sprintf('%s_%s',fbase,filename_coefs);
fpathf = fullfile(savedir,filename);
fprintf('Saving "%s" to "%s" ... ',filename,savedir); % must use newline symbol "\n" at the end
save(fpathf,'min_i1','min_sign1','boo1','xx1','Ind',...
    'in_Ghat','v_x1','Xu','nw','gu','n','covariate','vhats_n','vhats_p','minW',...
    'scaled_ATT_itr','ATE_itr','W_bs_Ghat_itr');

%% Output for table 
percent=sum(nw.*in_Ghat)/n;
savings=sum(gu.*in_Ghat)/n;
if (tcost)
    total_savings=savings*n*12;
else
    total_savings=savings*n*12/1000;
end
ci_lb=(minW-prctile(vhats_abs,95))/n;
ci_ub=(minW+prctile(vhats_abs,95))/n;

savings_EWM_onedim=nan(1,4);
savings_EWM_onedim=array2table(savings_EWM_onedim,'VariableNames',{'Covariates','Share\ treated','Net\ cost\ changes','Total\ cost\ changes'}); 
format short g
savings_EWM_onedim.(1) = {covariate};
savings_EWM_onedim.(2) =round(percent*100);
savings_EWM_onedim.(3) =round(savings,2);
savings_EWM_onedim.(4) =roundn(total_savings,3);

ci_lb=round(ci_lb,2);
ci_ub=round(ci_ub,2);

ci_table=nan(1,4);
ci_table=array2table(ci_table,'VariableNames',{'Covariates','Share\ treated','Net\ cost\ changes','Total\ cost\ changes'});
cilb_string=string(ci_lb);
ciub_string=string(ci_ub);
ci_for_table=strcat('(',cilb_string,',',ciub_string,')');
ci_table.(1)="";
ci_table.(2)="";
ci_table.(3)=ci_for_table;
ci_table.(4)="";
savings_EWM_onedim=[savings_EWM_onedim;ci_table];
end